#1.数据处理
data <- read.table("D:/geo-database-analysis/GSE44550/GPL17077-17467.txt",sep="\t",header=TRUE)
data <- read.table("D:/geo-database-analysis/GSE44550/GSE44550_series_matrix.txt",sep="\t",header=TRUE)
library(limma)
#kcol处输入gene symbol的所在列数
kcol=7
ann="D:/geo-database-analysis/GSE44550/GPL17077-17467.txt"
GSE="D:/geo-database-analysis/GSE44550/GSE44550_series_matrix.txt"
ann=read.table(ann,sep = "\t",header = F,check.names = F,fill = T)[,c(1,kcol)]
GSE=read.table(GSE,sep = "\t",header = T,check.names = F,fill = T)
af=merge(GSE,ann,by.x = colnames(GSE)[1],by.y = colnames(ann)[1])
rownames(af)=af[,ncol(af)]
af=af[,-c(1,ncol(af))]
af=as.data.frame(af)
af[,1:ncol(af)]=lapply(af[,1:ncol(af)],as.numeric)
af=as.matrix(af)
af=avereps(af)
af=af[order(af[,2]),]
rt=rbind(ID=colnames(af),af)
write.table(rt,file="matrix.txt",sep="\t",quote=F,col.names = F)
# 在 af=as.matrix(af) 之后，af=af[order(af[,2]),] 之前添加：
af <- af[rownames(af)!= "", ]  # 删除基因名为空的行

#画PCA图：看样本是否按“分组”聚类（比如癌vs正常组是否分开）
# 先转置矩阵（PCA要求行=样本，列=基因）
pca_data <- t(af)
# 做PCA
pca <- prcomp(pca_data, scale. = TRUE, center = TRUE)
# 提取前2个主成分（解释大部分差异）
pca_df <- as.data.frame(pca$x[, 1:2])
pca_df$Sample <- rownames(pca_df)
# 标注分组（示例：Group列，T=癌，N=正常）
pca_df$Group <- c("T", "N")  # 需根据实际样本类型修改！
# 画PCA图
library(ggplot2)
ggplot(pca_df, aes(x = PC1, y = PC2, color = Group, label = Sample)) +
  geom_point(size = 3) +
  geom_text(hjust = 1.2) +
  labs(title = "PCA of Samples", x = paste0("PC1 (", round(pca$sdev[1]^2/sum(pca$sdev^2)*100, 1), "%)"), 
       y = paste0("PC2 (", round(pca$sdev[2]^2/sum(pca$sdev^2)*100, 1), "%)")) +
  theme_bw()
#2.差异分析
# 加载必要的包
library(limma)
# 假设矩阵数据已经是af（行=基因，列=样本）
# 如果还没有加载数据，可以用以下代码读取matrix.txt（确保文件在工作目录）
# af <- read.table("matrix.txt", sep="\t", header=TRUE, row.names=1, check.names=FALSE)

# 查看矩阵结构（确认列是样本，行是基因）
head(af)
colnames(af)  # 应显示 "GSM1086291" "GSM1086292"

# 创建分组数据框（2个样本的分组）
group <- data.frame(
  sample = colnames(af),
  group = c("T", "N")  # 按实际分组修改，T=癌症，N=正常
)

# 构建设计矩阵（区分两组样本）
design <- model.matrix(~0 + factor(group$group))
colnames(design) <- c("N", "T")  # 列名对应分组

# 拟合线性模型
fit <- lmFit(af, design)

# 提取两组样本的表达量（假设第1列是T组，第2列是N组）
t_expr <- af[, 1]  # 癌症组表达量
n_expr <- af[, 2]  # 正常组表达量

# 计算log2倍数变化（T/N）
log2fc <- log2(t_expr / n_expr)

# 计算原始倍数变化（T/N）
fc <- t_expr / n_expr

# 构建结果数据框
diff_result <- data.frame(
  gene = rownames(af),
  T_expr = t_expr,       # 癌症组表达量
  N_expr = n_expr,       # 正常组表达量
  FC = fc,               # 倍数变化（T是N的多少倍）
  log2FC = log2fc        # log2转换的倍数变化
)

# 按log2FC绝对值排序（方便筛选差异大的基因）
diff_result <- diff_result[order(-abs(diff_result$log2FC)), ]

# 查看结果
head(diff_result, 10)

# 保存结果
write.table(diff_result, "diff_genes_fc.txt", sep="\t", quote=FALSE, row.names=FALSE)

# 筛选 log2FC > 1（即T是N的2倍以上，高表达）或 log2FC < -1（即N是T的2倍以上，低表达）的基因
significant_genes <- diff_result[abs(diff_result$log2FC) > 1, ]

# 查看显著差异基因
nrow(significant_genes)  # 数量
head(significant_genes)

#3.可视化分析（无 p 值）

#加载ggplot2包
library(ggplot2)
# 添加“是否显著”标签（基于FC阈值）
diff_result$significance <- ifelse(
  abs(diff_result$log2FC) > 1,
  ifelse(diff_result$log2FC > 1, "Up", "Down"),
  "Not significant"
)

#①绘制火山图（y轴用表达量均值，替代p值）
ggplot(diff_result, aes(x=log2FC, y=log2((T_expr + N_expr)/2), color=significance)) +
  geom_point(size=2) +
  scale_color_manual(values=c("blue", "red", "gray")) +
  geom_vline(xintercept=c(-1, 1), linetype="dashed") +  # 2倍差异阈值
  labs(title="Gene Expression Fold Change", 
       x="log2(Fold Change, T/N)", 
       y="log2(Average Expression)") +
  theme_bw()

#②绘制热图
library(pheatmap)

# 选择前30个差异最明显的基因（数量可调整）
top_genes <- head(diff_result$gene, 30)
top_af <- af[top_genes, ]

# 对基因表达量进行标准化（让不同基因的表达范围可比）
top_af_scaled <- t(scale(t(top_af)))  # 按行标准化

pheatmap(
  top_af_scaled,
  show_rownames = TRUE,  # 显示基因名
  show_colnames = TRUE,  # 显示样本名
  scale = "none",  # 已手动标准化，这里不重复
  color = colorRampPalette(c("green", "black", "red"))(100),  # 绿-黑-红表示低-中-高表达
  main = "Heatmap of Top Differentially Expressed Genes",
  annotation_col = data.frame(Group = c("T", "N")),  # 样本分组注释
  annotation_colors = list(Group = c("T" = "red", "N" = "blue"))  # 分组颜色
)

#③绘制散点图（x轴：N组，y轴：T组）

ggplot(diff_result, aes(x = N_expr, y = T_expr)) +
  geom_point(alpha = 0.6, color = ifelse(abs(diff_result$log2FC) > 1, "red", "gray50")) +
  geom_abline(slope = 1, intercept = 0, color = "blue", linetype = "dashed") +  # 对角线（表达量相等）
  scale_x_log10() +  # 对数转换，让分布更均匀
  scale_y_log10() +
  labs(title = "Gene Expression in T vs N", x = "Expression (N Group)", y = "Expression (T Group)") +
  theme_bw()

#④绘制MA图

# 计算每个基因的平均表达量（两个样本的均值）
diff_result$mean_expr <- (diff_result$T_expr + diff_result$N_expr) / 2
ggplot(diff_result, aes(x = log2(mean_expr), y = log2FC)) +
  geom_point(alpha = 0.7, color = ifelse(abs(diff_result$log2FC) > 1, "red", "gray50")) +
  geom_hline(yintercept = c(-1, 1), linetype = "dashed", color = "blue") +  # 2倍差异阈值
  labs(title = "MA Plot", x = "log2(Average Expression)", y = "log2(Fold Change)") +
  theme_bw()
